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1  Introduction 


Beginning  in  1984  and  extending  into  1993,  the  USAF  Armstrong  Laboratory, 
through  the  Aerospace  Medicine  Directorate’s  Clinical  Sciences  Division,  sup¬ 
ported  a  unique  effort  in  the  model-based  derivation  of  ECG  diagnostics;  the 
purpose  of  this  effort  was  development  of  an  enhanced  ECG  lead  set  and  di¬ 
agnostic  criteria  for  detection  of  asymptomatic  coronary  artery  disease  (CAD) 
in  USAF  aircrew  members.  This  work  was  principally  focused  in  a  contracted 
effort;  Dr.  Ronald  H.  Selvester  (University  of  Southern  California,  Rancho  Los 
Amigos  Campus,  and  Memorial  Heart  Institute,  Long  Beach  Memorial  Medical 
Center)  was  the  principal  investigator,  and  Col.  Gil  D.  Tolan  (retired)  was  the 
original  in-house  project  scientist.  Joseph  C.  Solomon,  the  co-investigator  and 
Dr.  Selvester’s  long-time  collaborator,  was  responsible  for  the  development  of 
the  simulation  software,  as  well  as  other  technical  matters,  as  required  by  the 
contract. 

Dr.  Selvester’s  work  on  the  development  of  ECG  diagnostic  criteria  has  been 
based  on  the  premise  that  an  advantage  could  be  gained  by  exploiting  knowledge 
derived  from  computer  simulation  of  the  underlsdng  cardiac  electrical  phenom¬ 
ena  ([1]  and  references  therein).  Validation  of  his  work  on  the  localization  and 
quantification  of  infarct,  using  criteria  developed  with  such  an  approach,  is  well 
documented  (references  in  [1],  [4],  and,  e.g.,  [5]).  The  effort  supported  by  the 
Clinical  Sciences  Division  was  intended  to  generalize  his  work  on  diagnostic 
criteria  for  infarct  to  the  case  of  ischemia. 

As  part  of  this  effort,  Selvester  and  Solomon  developed  a  computerized  for¬ 
ward  model  of  the  ECG;  it  consisted  of  a  ventricular  electrophysiology  (depo¬ 
larization  and  repolarization)  model,  a  multiple  fixed-dipole  equivalent  cardiac 
source  model,  a  torso  transfer  impedance  model,  and  a  model  for  derivation  of 
various  ECG  lead  systems  from  the  source  and  torso  models.  Details  of  Selvester 
and  Solomon’s  models  can  be  found  in  their  1985  Technical  Report  [1]. 

Due  to  funding  limitations,  their  models  were  not  developed  to  the  extent 
that  was  originally  intended,  and  did  not  incorporate  the  simulation  of  many 
aspects  of  normal  and  ischemic  electrophysiology  which  were  thought  to  be 
important  to  an  understanding  of  ECG  diagnostics  for  asymptomatic  CAD. 
Work  Unit  77552720  was  instituted  primarily  as  a  response  to  this  situation,  to 
provide  a  channel  through  which  in-house  work  could  be  performed  to  enhance 
the  contract  models  so  a  more  comprehensive  set  of  model-based  hypotheses 
could  be  generated.  The  two  principal  tasks  set  forth  at  the  inception  of  the 
work  unit  were  to: 

•  Generalize  the  simulation  model  of  ventricular  electrical  activity,  and 

•  Stabilize  the  numerical  behavior  of  the  torso  transfer  matrix  computation. 

These  two  issues  were  addressed  in  a  natural  progression.  Our  initial  work 
focused  on  development  of  an  enhanced  simulation  of  the  ventricular  depolar- 
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ization  and  repolarization  processes  which  would  have  the  capacity  to  simulate 
inhomogeneous  and  anisotropic  conduction,  re-entry  into  the  Purkinje  system, 
and  dispersion  of  refractoriness  in  a  way  that  could  be  easily  adapted  to  include 
other  electrophysiological  phenomena  which  might  arise  as  subjects  of  interest. 
The  resulting  simulation  is  phenomenological  rather  than  mechanistic,  requiring 
explicit  specification  of  action  potential  components  and  conduction  velocities 
throughout  the  ventricles;  on  the  other  hand,  one  has  complete  control  of  the 
simulation  scenario. 

As  an  outgrowth  of  an  attempt  to  validate  the  ventricular  model  in  terms 
of  internal  consistency  with  the  entire  ECG  forward  model  scheme,  we  began 
work  on  the  problem  of  stabilizing  the  torso  transfer  matrix  computation,  but 
this  work  was  not  completed. 

In  this  technical  paper,  we  discuss  only  our  discrete-event  simulation  model 
of  myocardial  electrical  activity.  It  is  the  only  complete  component  of  the  work 
intended  to  be  done  under  the  work  unit,  and  its  application  provided  the 
impetus  for  other  work  that  was  attempted  but  not  completed.  Following  a 
description  of  the  simulation  model,  we  describe  the  general  results  and  conclude 
with  a  discussion  of  how  this  work  relates  to  broader  issues,  as  well  as  “lessons 
learned”  from  the  effort. 


3 


2  The  Discrete-Event  Simulation  Model 

Like  Selvester  and  Solomon’s  model  of  the  heart,  our  model  is  a  “simulation 
model”  rather  than  as  a  “mathematical  model”  because  it  is  algorithmic  in 
nature,  as  opposed  to  mathematical.  In  other  words,  it  is  defined  by  the  al¬ 
gorithms  embodied  in  the  simulation  software  rather  than  by  a  set  of  partial 
differential  equations,  as  would  be  the  case  if  we  were  modeling  the  ion  kinetics. 
Likewise,  it  is  a  simulation  of  bioelectric  phenomena  rather  than  a  model  of 
electrophysiological  mechanism. 

The  core  of  Selvester  and  Solomon’s  model,  as  well  as  ours,  is  the  algo¬ 
rithm  for  computing  the  depolarization  wavefront.  Conceptually,  Selvester  and 
Solomon  base  their  algorithm  on  Huygens’  method  of  wavefront  propagation  [1]. 
In  this  method,  the  wavefront  at  time  t2  is  formed  from  the  wavefront  at  a  pre¬ 
vious  elementary  time  step,  ti,  by  forming  the  envelope  of  elementary  wavelets 
emanating  from  each  point  of  the  waveficont  at  tune  ti,  where  the  elementary 
wavelets  depend  on  a  uniform,  direction-independent  conduction  velocity.  Our 
model  is  conceptually  based  on  a  generalization  of  Huygens’  method  in  which 
each  point  of  the  underlying  space  (e.g.,  the  ventricular  myocardium)  possesses 
an  “indicatrbc”  of  direction-dependent  conduction  velocities  [6].  In  both  models, 
the  cells  are  conceived  of  as  being  uncoupled  during  repolarization. 

A  major  goal  of  our  work  was  to  maintain  complete  input-output  data  com¬ 
patibility  with  the  Selvester-Solomon  simulation,  while  generalizing  the  class 
of  electrophysiological  phenomena  that  could  be  simulated.  Data  compatibility 
was  achieved,  but  generalization  of  the  simulation  required  a  complete  redesign 
of  the  software. 

Acknowledgement  of  Data  Source.  All  of  the  input  data  for  our  model  was 
provided  under  contract  along  with  the  Selvester-Solomon  simulation  software. 
AU  of  the  model  elements  described  below  were  implemented  using  this  data, 
which  was  developed  by  Joe  Solomon  under  the  guidance  of  Dr.  Selvester. 
Additional  details  on  many  of  these  data  elements  can  be  found  in  the  technical 
report  by  Selvester  and  Solomon  [1]. 

2.1  Anatomical  and  Electrophysiological  Structure 

The  global,  or  large-scale,  anatomical  structure  is  defined  by  a  digital  represen¬ 
tation  of  an  actual  ventricular  geometry,  at  diastole,  obtained  under  contract 
from  Selvester  and  Solomon  [1].  Global  coordinates  are  assigned  to  the  points 
of  the  geometry  by  setting  it  in  the  positive  octant  of  the  conventional  three- 
dimensional,  orthogonal,  right-handed  x-y-z  coordinate  system  with  axes  scaled 
to  measure  in  1mm  units,  and  with  the  z-axis  measuring  the  extent  from  apex  to 
base.  So  defined,  the  geometry  fits  easily  within  a  128mm3  cube.  In  its  priirutive 
form,  this  geometry  is  defined  by  occupancy  of  points  on  the  three-dimensional 
integer  lattice  (  points  (i,  j,  k),  where  i,  j,  and  k  are  integers)  within  this  cube. 


4 


Conceptually,  we  take  these  occupied  points  to  represent  the  center  points  of  1 
mm^  blocks  of  cardiac  tissue  which  make  up  the  elementary  cellular  units  of  the 
model. 

The  ventricular  structure  is  decomposed  into  two  basic  types  of  spatially 
discrete  cellular  units:  Purkinje  units  and  myocardial  units.  Purkinje  units 
represent  blocks  of  tissue  with  a  dense  Purkinje  fiber  distribution  and  myocardial 
units  represent  blocks  predominantly  incorporating  myocardium.  Each  cellular 
unit  is  associated  with  a  set  of  electrophysiological  parameters  and  properties 
which  are  used  in  the  simulation  to  reproduce  desired  bioelectric  phenomena. 
The  geometric  discretization  embodies  an  implicit  assumption  that  the  small 
anatomical  region  represented  by  each  cellular  unit  is  homogeneous  enough  that 
the  corresponding  parameters  and  properties  of  its  biological  constituents  does 
not  vary  significantly  from  the  “average  values”  assigned  to  the  cellular  unit. 

The  ventricrdar  tissue  model  is  represented  in  the  simulation  model  by  two 
data  structures: 

•  A  128  X  128  X  128  array,  each  element  of  which  contains  a  pointer.  If 
the  element  corresponds  to  coordinates  for  a  ventricular  cellular  unit,  the 
pointer  indicates  an  element  in  an  array  of  data  structures  containing  rep¬ 
resentations  for  the  electrophysiological  components  of  the  cellular  units. 

If  the  element  does  not  correspond  to  coordinates  for  a  ventricular  cellular 
unit,  the  pointer  indicates  a  “null”  value. 

•  An  array  of  data  structures  containing  the  specifications  for  the  cellu¬ 
lar  units  in  terms  of  electrophysiological  components,  such  as  conduction 
properties  and  action  potenticil  parameters. 

The  rationale  for  this  dual  structure  is  that  the  full  ventricular  geometry  con¬ 
tains  287, 140  myocardial  units  and  56, 090  Purkinje  units,  for  a  total  of  343,  230 
cellular  units,  comprising  approximately  16.4%  of  the  128®  coordinate  array.  As¬ 
signment  of  a  data  structure  to  each  point  of  the  array  would  have  resulted  in 
an  enormous,  unnecessary  consumption  of  computer  memory.  In  fact,  an  earlier 
design  incorporating  such  an  assignment  was  tested  on  the  Digital  Equipment 
Corporation  (DEC)  VAX  computers  available  for  the  simulations.  The  intense 
referencing  of  this  aggregate  data  structure  by  our  simulation  algorithm  pro¬ 
duced  severe  memory  paging  problems  in  the  VAX  computer,  resulting  in  an 
extremely  impractical  program  runtime. 

The  data  structure  for  each  cellular  unit  contains  values  for  the  following 
electrophysiological  components: 

•  Cell  Type  -  Purkinje  or  myocardium. 

•  Conduction  Velocity  Indicatrix  -  specification  of  possibly  direction-dependent 
conduction  velocity,  relative  to  the  center-location  of  the  cellular  unit,  in 

a  way  that  is  commensurate  with  the  depolarization  wave  propagation 
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algorithm  (see  below).  The  specification  in  terms  of  conduction  velocity 
is  replaceable  by  liber  orientation  and  conductivity  tensor,  depending  on 
the  algorithm  for  computing  propagation. 

•  Action  Potential  Waveform  -  either  a  pointer  to  a  regional  action  poten¬ 
tial  parameter  set  or  an  individual  set  of  action  potential  parameters, 
depending  on  the  configuration,  or  objective,  of  the  simulation.  The  ac¬ 
tion  potential  function  is  a  continuous-time  version  of  an  ad  hoc  formula 
derived  by  Joe  Solomon  [1].  Its  functional  form  is  fixed  for  all  cells,  but 
its  morphology  can  be  varied  by  modulation  of  a  number  of  parameters. 
The  (default)  distributions  of  these  parameters,  across  the  cellular  units, 
are  specified  to  approximate  the  spatial  distribution  of  action  potential 
morphologies  in  a  normal  human  heart.  These  parameters  include 

—  Resting  potential, 

—  Rise-time  of  action  potential, 

—  Depolarization  potential, 

—  Absolute  refractory  period,  and 
—  Recovery  rate. 

•  Action  Potential  Firing  Time  -  the  time,  in  milliseconds(ms),  at  which  the 
cellular  unit  was  last  activated.  Time  is  typically  measured  in  reference  to 
the  first  activation  of  a  Purkinje  or  myocardial  unit  being  designated  as 
time  zero  (0).  The  simulation  begins  with  the  action  potential  firing  time 
of  all  cellular  units  set  to  a  common  negative  value  such  that  they  are  in 
a  normal  resting  state  at  time  zero. 

Note.  In  the  following,  for  the  sake  of  simphcity,  we  use  the  term  “cell”  in¬ 
stead  of  “cellular  unit”  to  represent  the  macroscopic  block  of  tissue  previously 
described. 

2.2  Propagation  of  the  Depolarization  Wavefront  and  Re¬ 
polarization 

Briefly,  the  simulation  model  begins  with  initial  activation  of  bundle  entry  sites 
among  the  cells  designated  eis  Purkinje  type.  From  this  beginning,  activation 
is  propagated  by  electrotonic  interaction  between  cells,  as  well  as  by  additional 
bundle  entries,  until  the  entire  ventricular  mass  is  activated.  Activation  of  a  cell 
marks  the  temporal  starting  point  for  its  associated  action  potential  waveform, 
from  which  the  course  of  repolarization  is  computed.  The  complexity  in  this 
process  arises  from  computation  of  the  depolarization  wavefront. 

The  initial  activation  of  bundle  entry  sites  results  in  the  formation  of  an 
initial  depolarization  wavefront.  Activation  times  are  computed  for  aU  repolar¬ 
ized  cellular  units  within  the  26  cell  neighborhood  of  each  of  the  cells  in  the 
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wavefront.  The  algorithm  for  computation  of  the  activation  times  accounts  for 
the  effects  of  a  number  of  anatomical  and  electrophysiological  phenomena,  as 
described  below.  During  the  computations  for  each  wavefront,  the  neighboring 
cells  are  “scheduled”  and  “rescheduled”  to  insure  that  their  activation  times  are 
the  earhest  times  at  which  propagation  would  occur  relative  to  the  wavefront. 
The  next  wavefront  consists  of  all  those  cells  having  the  “next”  activation  time 
appearing  in  the  activation-event  scheduling  structure.  In  this  process,  there  is 
no  pre-specified  minimal  time  step.  The  interplay  of  Purkinje  and  myocardial 
conduction  velocities  with  three-dimensional  geometric  relationships  results  in 
a  dense  structure  of  activation  times  which  have  to  be  managed  with  special 
data  structures. 

Our  simulation  model  has  the  potential  of  reproducing  the  following  phe¬ 
nomenological  features  of  myocardial  activation: 

•  Electrotonic  Interaction  Between  Cellular  Units. 

Propagation  of  the  depolarization  wavefront  proceeds  from  ceU  to  cell, 
with  an  activated  ceU  propagating  the  wavefront  to  a  neighboring  ceU  at 
the  end  of  a  time  interval  depending  on  the  conduction  velocity  at  the 
activated  ceU  and  the  distance  to  the  neighboring  ceU. 

The  neighborhood  of  a  given  ceU  consists  of  all  ceUs  whose  centers  are 
strictly  less  than  2  mm  from  the  center  of  the  given  ceU.  Defined  in  this 
way,  the  neighborhood  of  a  given  ceU  can  be  visualized  as  containing  aU 
ceUs  whose  centers  lie  on  a  face,  on  an  edge,  or  at  a  vertex  of  the  cube 
centered  at  the  given  ceU  and  having  edges  of  2  mm  in  length.  Thus,  the 
neighborhood  of  a  ceU  contains  at  most  26  other  cells,  depending  on  how 
many  of  the  surrounding  lattice  points  are  occupied  by  ceUular  units  and 
how  many  surrounding  lattice  points  correspond  to  the  space  outside  of 
the  heart. 

In  the  case  of  a  homogeneous,  isotropic  conduction  velocity,  of  magnitude 
V,  within  the  region  of  computation,  the  propagation  time  from  an  acti¬ 
vated  ceU  located  at  coordinate  vector  ®  to  a  neighboring  ceU  located  at 
coordinate  vector  y  is  given  by  the  distance  between  the  ceUs  divided  by 
the  magnitude  of  the  velocity: 

T(®,  =  \y-  xl/v. 

•  Inhomogeneity  in  Cell  Type. 

There  are  two  types  of  cells  in  our  model,  resulting  in  three  different  possi¬ 
ble  interactions:  Purkinje-Purkinje,  Purkinje-myocardium,  and  myocardium- 
myocardium.  In  principle,  propagation  between  ceUs  of  different  type 
should  be  modeled  differently  than  propagation  among  ceUs  of  the  same 
type.  However,  in  our  simulation,  the  Purkinje  ceUular  units  are  taken  to 
model  an  averaged  behavior  for  myocardium  which  is  densely  interlaced 
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with  Purkinje  fibers.  Therefore,  a  Purkinje-to-myocardium  interaction 
(and  vice  versa)  is  defined  as  equivalent  to  a  myocardium-to-myocardium 
interaction.  On  the  other  hand,  we  foUow  Selvester  and  Solomon  in  spec¬ 
ifying  that  the  propagation  in  a  Purkinje-Purkinje  interaction  be  three 
times  faster  than  among  myocardial  cells.  (This  ratio  varies  significantly 
among  models  reported  in  the  research  literature.)  Consequently,  we 
model  Purkinje  cells  with  a  conduction  velocity  indicatrix  of  the  same 
relative  magnitude  as  myocardium,  but  increase  its  magnitude  by  a  factor 
of  three  when  applying  it  to  a  neighboring  Purkinje  cell. 

•  Inhomogeneity  in  Conduction  Velocity. 

Deviations  in  conduction  velocity  are  known  to  occur  with  ischemia;  thus, 
an  essential  feature  of  the  model  is  the  ability  to  account  for  inhomo¬ 
geneities  in  conduction  velocity  within  the  myocardium. 

The  structure  of  our  tissue  model  allows  us  to  specify  a  different  con¬ 
duction  velocity  indicatrix  for  each  cell  in  the  ventricular  anatomy.  Our 
simulations  took  the  form  of  either  scaling  the  normal  myocardial  conduc¬ 
tion  velocities  in  a  region  by  a  fixed  factor  (e.g.,  one-half)  or  randomly 
modulating  the  normal  conduction  velocities  in  a  region. 

•  Fiber  Anatomy  and  Anisotropic  Conduction. 

Individual  myocardial  cells  are  elongated  and  connected  in  a  roughly  end- 
to-end  manner  to  form  fibers,  which  in  turn  serve  as  components  of  larger 
structures.  These  structures  are  organized  into  a  complex  architecture 
which  is  optimized,  especially  in  the  left  ventricle,  for  efficient  contraction 
and  emptying  of  the  ventricular  chambers. 

The  fiber  structures  of  the  heart  also  play  a  significant  roles  in  the  forma¬ 
tion  of  electrophysiological  phenomena.  The  special  nature  of  the  junc¬ 
tions  between  cells  in  a  fiber  results  in  higher  electrical  conductivity  along 
the  length  of  the  fiber  than  across  it.  As  a  result,  the  conduction  velocity 
of  the  depolarization  wavefront  is  directionally  dependent,  or  anisotropic, 
depending  on  the  local  average  fiber  direction. 

At  the  global  level,  the  structure  remains  complex,  but  generally  results 
in  a  circumferential  conduction  that  is  faster  than  the  transmural  con¬ 
duction.  It  is  thought  by  many  researchers  that  the  anisotropic  nature 
of  the  conduction  does  not  affect  the  normal  activation  sequence  for  the 
ventricles,  because  the  propagation  of  the  depolarization  wavefront  gen¬ 
erally  takes  a  transmural  course  from  endocardium  to  epicardium.  On 
the  other  hand,  the  anisotropic  nature  of  the  conduction  velocities  would 
seem  important  in  accurately  modeling  cases  in  which  the  depolarization 
wavefront  does  not  propagate  transmurally  in  a  uniform  manner,  such  as 
in  some  cases  of  infarct  and  ischemia,  bundle  branch  blocks,  and  many 
arrhythmias. 
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Because  of  the  simulation  model’s  organization,  the  mechanism  of  conduc¬ 
tion  can  be  specified  by  the  user  in  a  manner  most  convenient  for  the  pur¬ 
pose  at  hand.  In  most  published  simulation  studies  involving  anisotropic 
conduction  ([7]  and  references  therein),  actual  or  approximate  fiber  di¬ 
rections  are  used  to  define  the  directional  dependence  of  electrical  con¬ 
duction,  which  is  in  turn  used  to  define  conduction  velocities  in  some  ad 
hoc  manner.  We  took  a  simpler  approach  and  worked  directly  in  terms  of 
conduction  velocities.  The  natural  way  to  define  a  fiber- based  conduction 
velocity  is  to  specify  a  conduction  velocity  along  the  length  of  the  fiber 
and  a  circumferential  velocity  along  any  orthogonal  axis.  From  this,  a 
fiber-based  conduction  velocity  can  be  computed  for  any  direction.  On 
the  other  hand,  we  had  already  defined  a  global  coordinate  system  for 
specifying  the  location  of  cells,  so  we  preferred  to  define  our  conduction 
velocities  directly  in  terms  of  the  same  global  coordinate  cixis  directions. 

Specifically,  the  conduction  velocity  indicatrix  for  a  given  cell  at  coordinate 
vector  X  can  be  specified  by  a  velocity  vector,  r,  which  might  be  formed 
by  projecting  the  longitudinal  and  circumferential  components  of  a  fiber- 
based  conduction  velocity  onto  rectangular  axes  originating  at  x.  Then, 
the  propagation  time  from  the  given  cell  to  a  neighboring  ceU  centered  at  y 
is  defined  as  the  distance  between  the  two  cells  divided  by  the  magnitude 
of  the  component  of  the  conduction  velocity  projected  onto  the  direction 
indicator  from  x  to  y-. 


T{y,  x)  = 


\y-  s| 


V  o 


\y-^\ 


) 


where  “o”  denotes  an  inner  product  of  vectors. 


•  State-Dependent  Propagation  between  Cells. 

The  excitability  (susceptibility  to  depolarization)  of  a  myocardial  cell  is 
dependent  on  its  “refractory  state,”  which  can  be  thought  of  as  a  state  cor¬ 
responding  to  certain  conditions  at  the  ionic  level.  In  the  fully  repolarized 
or  resting  state,  the  ceU  is  capable  of  participating  in  normal  ionic  current 
flow,  with  the  capability  for  normal  conduction  and  production  of  a  nor¬ 
mal  action  potential.  CeU  depolarization  and  recovery  foUows  a  sequence 
of  ionic  current  flows  which  are  interuptable  only  under  certain  conditions. 
When  the  action  potential  of  a  myocardial  ceU  is  in  the  plateau  phase  of 
its  waveform,  it  is  in  an  “absolute  refractory  state”  and  is  not  interupt¬ 
able.  When  it  is  repolarizing  in  the  downward  slope  of  the  waveform, 
however,  the  underlying  ionic  current  flows  are  partiaUy  interuptable,  and 
depolarization  can  again  be  achieved.  In  this  “relative  refractory  state”, 
the  response  is  muted,  with  decreased  conduction  velocity  and  distortion 
of  the  resulting  action  potential. 
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In  our  model,  this  phenomenon  is  simulated  by  modulating  the  conduc¬ 
tion  velocity  between  an  activated  cell  and  a  neighboring  cell  according 
to  the  refractory  state  of  the  neighboring  cell.  If  the  neighboring  cell 
is  completely  repolarized,  then  the  propagation  time  between  the  cells  is 
computed  using  the  normal  conduction  velocity  indicatrix  of  the  activated 
cell.  If  the  neighboring  cell  is  in  an  absolute  refractory  state  (i.e.,  within 
the  absolute  refractory  period  as  specified  for  its  action  potential  wave¬ 
form)  it  cannot  be  reactivated,  so  the  conduction  velocity  between  the 
cells  is  set  to  zero.  If  the  neighboring  cell  is  in  a  relative  refractory  state 
(i.e.,  beyond  its  absolute  refractory  period,  but  not  yet  fuUy  repolarized) 
then  the  conduction  velocity  indicatrix  of  the  activated  cell  is  scaled  by 
the  following  ratio: 

Action  Potential  —  Depolarization  Potential 
Resting  Potential  —  Depolarization  Potential 

We  used  this  scheme  primarily  for  proof  of  concept.  More  realistic  func¬ 
tional  relationships  can  be  found  by  consulting  the  research  literature  on 
refractory  states  of  excitable  cells. 

In  our  implementations,  we  only  modulated  the  conduction  velocities,  and 
did  not  attempt  to  modulate  the  action  potential  morphologies,  although 
its  incorporation  into  the  software  should  be  easily  managed  for  a  model 
specified  in  terms  of  functional  relationships. 

This  is  not  the  whole  story.  Of  particular  importance  in  modeling  arrhythmias, 
as  well  as  exercise  effects  (there  are  more  difficult  problems  in  this  case),  is  the 
dependency  of  conduction  and  action  potential  morphology  on  the  frequency  of 
activation.  Genereil  functional  forms  for  some  of  these  dependencies  are  available 
from  the  research  literature,  but  we  did  not  explore  such  a  capability  as  part  of 
our  work. 

The  Simulation  Algorithm.  The  main  data  components  and  behavioral  el¬ 
ements  of  the  simulation  model  have  been  described  above.  We  will  now  present 
a  succinct  description  of  the  most  general  simulation  algorithm  we  implemented. 
Because  of  the  admitted  provisional  nature  of  some  of  the  modeling  components, 
this  should  be  considered  more  of  an  example  than  a  definitive  exposition. 

Let  VM  denote  the  ventricular  geometry,  including  both  Purkinje  and  my¬ 
ocardial  cells,  and  let  WF  denote  the  depolarization  wavefront.  SCHEDULE 
is  a  data  structure  which  is  used  to  accumulate  the  “next”  depolarization  wave- 
front  by  managing  “predicted”  firing  (depolarization)  times  for  cells  neighboring 
the  wavefronts,  and  updating  them  as  necessary.  The  simulation  is  initialized 
by  constructing  the  ventricular  geometry  and  array  of  cell  parameters,  and  by 
initializing  SCHEDULE  with  a  set  of  locations  and  firing  times  for  entry  sites 
of  the  conduction  bundles  into  the  Purkinje  system.  In  the  following  example. 
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we  shall  reference  a  cell  in  the  ventricular  structure  by  its  location  coordinate 
vector  (e.g.,  x). 

The  basic  algorithm  is  as  follows; 

•  Initialize  the  simulation. 

•  Do  the  following,  until  the  myocardium  is  fully  repolarized: 

—  Extract  the  next  wave&ont,  WF,  at  time  T,  from  SCHEDULE. 


—  Propagate  the  wavefront,  WF: 

For  every  £  in  WF,  compute  the  local  wavelet  propagation: 


*  Find  the  neighborhood  of  x,  N{x)  in  FAf. 

*  Determine  the  tissue  type  of  x,  TYPE{x). 


*  for  every  j/ in  N{x): 

•  Determine  the  tissue  type  of  y,  TYPE{‘^,  and  compute 
the  type  factor  TF,  where  TF  =  Z  if  both  TYPE[x)  and 
TYPE{y)  are  Purkinje,  and  TF  =  1  otherwise. 


•  Determine  the  state  of  y: 


STATE{y)  = 


AP{y,T)-Va 

V,-Vi 


where  AP{y,T)  is  the  action  potential  voltage  of  cell  y  at 
time  T  measured  relative  to  the  last  action  potential  firing 
time  for  y ,  is  the  depolarization  potential  of  y,  and  Vr  is 
the  resting  potential  of  y. 


Compute  the  conduction  velocity  at  x  in  the  direction  of  y, 
as  dependent  on  cell  type  and  refractory  state: 

V{x,^  ^TFx  STATE{f)  x  . 


-  Compute  the  predicted  firing  time,  FT ,  of  y: 


FT  =  T  + 


|y-»l 

V{x,^' 
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•  Schedule  y  to  fire  at  time  FT: 

PUSEJSCHEDULE{y,  FT). 

Note  that  the  entire  section  following  “For  every  ®  in  WF:”  is  encapsulated  in 
a  single  subroutine  of  the  simulation  software,  acting  on  the  tissue  model,  and 
could  easily  be  replaced  by  more  general  functions. 


2.3  Software  Engineering 

There  are  two  computational  aspects  of  our  simulation  software  which  are  worth 
noting,  although  their  discussion  is  not  essential  for  understanding  of  the  con¬ 
ceptual  model  itself. 

First,  the  simulation  model  was  initially  regarded  as  being  provisional.  We 
expected  to  upgrade  various  components  as  time  passed,  so  we  constructed  the 
software  to  facilitate  such  occurrences  by  modularizing  at  the  level  where  we 
expected  the  changes  to  take  place.  This  led  us  to  consider  the  simulation  im¬ 
plementation  in  terms  of  correspondences  between  the  conceptual  components 
of  the  model  and  data  objects  in  the  software;  in  other  words,  form  the  view¬ 
point  of  an  object-oriented  design  philosophy  (see  [8]  and  references  therein). 
On  the  other  hand,  to  ensure  software  compatibility  with  the  contracted  work 
of  Selvester  and  Solomon,  we  were  implementing  our  simulation  in  Fortran  77 , 
which  does  not  support  object-oriented  programming  techniques.  Therefore,  we 
adopted  Fortran  programming  conventions,  some  of  which  were  drawn  ftom  the 
literature  ([9],  [10],  [8]),  to  implement  a  weak  form  of  data  abstraction  (encap¬ 
sulation  of  data  and  the  subroutines  that  operate  upon  it).  This  approach  was 
adequate  for  our  needs,  and  considered  a  success;  it  enabled  us  to  make  sev¬ 
eral  major  changes  to  specific  software  components  during  development  without 
requiring  significant  changes  to  other  components.  In  terms  of  the  conceptual 
components  of  the  model,  another  advantage  was  the  resulting  structural  co¬ 
herence  of  the  software,  which  made  the  code  easier  to  understand. 

The  second  computational  point  of  interest  is  in  our  choice  of  event-scheduling 
methodology.  While  there  have  been  a  number  of  reports  in  the  literature  of 
simulations  based  on  the  methodology  of  Selvester  and  Solomon’s  early  work 
(recently:  [7],  [11]),  few  simulations  of  cardiac  electrophysiology  have  been  im¬ 
plemented  using  a  discrete-event  simulation  approach.  The  discrete-event  simu¬ 
lations  that  have  been  reported  are  usually  implemented  using  specialized  data 
structures  and  algorithms  for  managing  the  timing  of  cell  activations  ([12],  [13]). 
We  wanted  to  avoid  such  an  approach  because  it  seemed  to  inhibit  flexibility  for 
future  software  evolution  and  would  place  the  work  outside  of  the  mainstream 
of  simulation  technology.  Such  a  decision,  however,  left  us  needing  methods  to 
efficiently  handle  the  scheduling  of  many  thousands  of  events. 

Using  standard  linked-lists  was  untenably  slow,  so  we  reviewed  the  literature 
and  found  a  study  by  Jones  which  reported  empirical  comparisons  of  priority- 
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queues  [14].  The  best  performance  for  very  large  queues,  under  most  of  the 
test  models,  was  obtained  by  the  splay  trees  of  Sleator  and  Tarjan  [15].  The 
implementation  of  splay  tree  scheduling  algorithms  resulted  in  much  faster  sim¬ 
ulations,  which  could  be  completed  overnight  when  run  in  batch  mode  on  the 
available  DEC  VAX  computers.  The  details  of  the  basic  algorithm  are  complex, 
and  would  best  be  attained  from  the  paper  of  Sleator  and  Tarjan,  where  they 
are  presented  in  a  “pseudocode”  format.  We  only  implemented  the  most  basic 
version,  although  there  are  a  number  of  possible  optimizations  cited  in  the  their 
paper. 

Our  initial  implementation  used  splay  trees  for  managing  events  using  two 
search  keys;  first  events  are  ordered  according  to  a  time  key,  and  then  for  each 
time  they  are  ordered  according  to  a  coded  representation  of  the  cell’s  spatial 
coordinates  in  the  128^  lattice.  This  algorithm  was  elegant  in  structure,  but 
required  additional  computations  for  the  coding  and  decoding  of  coordinates. 
Our  final  method  was  a  hybrid  approach,  using  a  splay  tree  for  management 
of  the  time  key,  and  time-indexed  bins,  implemented  as  linked  lists,  to  index 
the  cells  sched^ed  to  fire  at  each  time  key.  This  latter  component  was  derived 
from  a  suggestion  offered  independently  by  Col.  Gil  Tolan  and  Dr.  Sherwood 
Samn,  and  represents  a  successful  trade-off  of  computer  memory  for  speed  of 
computation. 
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3  Results 


The  simulation  model  was  first  tested  by  running  simulations  for  two-dimensional 
arrays  of  the  myocardial  cellular  units  to  confirm  that  the  algorithmic  machin¬ 
ery  for  computing  inhomogeneous  and  anisotropic  conduction  would  behave  as 
intended.  A  number  of  simple  conduction  scenarios  were  fabricated  to  exer¬ 
cise  the  algorithm.  The  results  were  displayed  graphically  on  a  Tektronix  4111 
computer  graphics  terminal. 

Next,  we  tested  the  fuU  three-dimensional  ventricular  geometry  model.  The 
output  from  these  simulations  was  either  an  interactive  computer  graphics  dis¬ 
play  or  a  batch-computed  animation  sequence  for  display  on  a  Siemens  DELTAVi- 
sion  medical  imaging  system.  The  display  format  was  basically  an  enhanced 
version  of  the  design  developed  under  contract  by  Selvester  and  Solomon  (see 
[16]  for  an  example).  Using  these  displays,  the  progression  of  depolarization 
and  repolarization  of  the  ventricles  could  be  directly  observed.  For  this  type  of 
simulation,  in  which  the  results  of  the  computation  are  distributed  in  space  as 
well  as  in  time,  computer  graphics  representation  of  the  results  seems  to  be  one 
of  the  few  direct  ways  to  confirm  correct  behavior  of  the  model.  We  are  unable 
to  include  images  from  the  computer  graphics  displays  in  this  report. 

Our  primary  design  constraint  in  developing  this  simulation  was  compatibil¬ 
ity  with  the  modeling  framework  of  Selvester  and  Solomon.  Our  first  test  using 
the  ventricular  geometry  was  duplication  of  the  simulations  of  their  ventricular 
model  (homogeneous,  isotropic  conduction  throughout  the  myocardium,  with 
conduction  in  the  Purkinje  system  specified  as  three  times  faster),  for  normal 
and  infarct  cases.  This  was  accomplished  by  setting  our  model’s  conduction 
parameters  at  equivalent  veilues  to  those  implicitly  found  in  their  model.  We 
found  excellent  agreement  between  the  two  models,  as  expected,  since  Selvester 
and  Solomon’s  model  served  as  the  default  behavior  of  our  model. 

To  test  the  mechanism  for  computing  inhomogeneous  conduction  in  its  in¬ 
tended  setting,  we  used  the  infarct  geometry  from  Selvester  and  Solomon’s 
model  to  define  a  region  of  slowed  conduction  (as  in  a  hypothetical  ischemic  re¬ 
gion)  in  the  ventricular  myocardium.  For  test  purposes,  the  slower  conduction 
velocity  was  simply  one-half  normal  velocity.  The  computer  graphics  display 
revealed  the  expected  lag  in  the  depolarization  wavefront  as  it  passed  through 
the  affected  region. 

Finally,  to  further  test  inhomogeneous  conduction,  as  well  as  propagation 
between  different  cell  types,  we  simulated  a  full  right  bundle  branch  block  with 
activation  originating  in  the  left  bundle,  as  well  as  a  similar  simulation  for  the 
case  of  full  left  bundle  branch  block.  The  depolarization  wavefront  proceeded 
as  expected,  passing  across  the  septum  and  re-entering  the  Purkinje  network 
below  the  block. 

These  scenarios  are  simplistic,  but  are  merely  intended  to  be  tests  and 
demonstrations  of  the  simulation  software  capabilities.  The  realism  of  the  simu¬ 
lations  would  be  expected  to  increase  with  increasing  realism  in  the  specification 
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of  input  parameters.  However,  we  were  unable  to  attempt  to  refine  these  pa¬ 
rameters. 
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4  Discussions 

4.1  The  Simulation  Model  in  a  Broader  Context 

Early  in  the  development  of  our  simulation  model,  we  began  to  view  it  as  more 
of  a  framework  for  model  implementation  than  as  a  single  fixed  model.  We  ex¬ 
pected  to  change  the  model’s  inner  mechanisms  as  it  was  applied  to  the  analysis 
of  increasingly  detailed  questions  regarding  the  electro  cardiology  of  ischemia. 
In  this  sense,  the  rudimentary  form  of  data  abstraction  which  we  applied  in 
development  of  our  software  was  important.  It  explicitly  organized  the  software 
into  a  framework  where  components  could  be  changed  with  little  impact  on  the 
rest  of  the  program. 

In  light  of  its  extendibUity,  we  consider  this  simulation  framework  to  be  as 
capable  as  any  phenomenological  simulation  model  of  the  ventricles  reported  in 
the  research  literature.  As  a  specific  model  of  ventricular  electrical  activity,  it 
lacks  a  complete  specification  of  realistic  electrophysiological  parameters,  as  well 
as  complete  development  of  certain  functional  characteristics  of  propagation, 
such  as  Jiequency  dependence  of  conduction  velocity. 

We  also  preferred  to  view  this  simulation  as  a  framework  rather  than  a 
fixed  model  because  the  simulation  of  “abnormal”  electrophysiology  in  terms  of 
macroscopic  phenomena  is  fundamentally  flawed.  Our  development  of  a  discrete, 
phenomenological  model  was  initially  mandated  by  a  desire  for  compatibility 
with  Selvester  and  Solomon’s  work,  as  well  as  by  the  fact  that  it  was  simply 
what  we  understood  how  to  do  at  the  time.  As  our  perspective  evolved  and  our 
knowledge  grew,  we  began  to  more  fully  appreciate  the  advantages  of  directly 
modeling  the  biophysical  mechanisms  underlying  the  phenomena.  In  fact,  an 
approach  to  the  problem  on  the  biophysical  level  is  essential  if  one  wishes  to 
begin  to  understand  what  is  happening  during  ischemia,  during  exercise,  or 
during  interactions  with  mechanical  fields  (e.g.,  myocardial  contraction  with 
and  without  -fGz  stress).  On  the  other  hand,  simulation  in  terms  of  biophysics 
is  not  tenable  at  the  scde  of  the  ventricular  muscle  mass.  In  fact,  simulations 
using  detailed  ionic  current  models  at  the  scale  of  the  whole  heart  are  beyond 
the  capabilities  of  current  supercomputers  for  known  methods  of  integration. 

Consequently,  strategies  which  balance  the  need  for  electrophysiological  de¬ 
tail  against  the  computational  burden  of  simulating  phenomena  at  the  scale 
of  the  ventricular  muscle  mass  are  needed.  Discrete,  phenomenological  simu¬ 
lation  may  yet  be  useful  in  the  development  of  such  strategies.  Specifically, 
studies  involving  detailed  biophysical  models  can  be  used  to  derive  data  and 
phenomenologiccil  models  of  dynamics  for  use  in  discrete-event  simulations  on 
the  large  scale.  After  realizing  this  at  a  fairly  early  stage,  and  beginning  work  in 
this  context,  we  were  pleased  to  find  these  sentiments  reflected  in  comments  in 
the  recent  paper  of  Manor,  et  al.,  on  propagation  of  neuronal  action  potentials 
[17]: 

...  it  is  clear  that  the  [detailed]  modeling  approach  is  limited  when 
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one  wishes  to  model  large  ...  systems  for  long  periods  of  time  ...  In 
such  cases  some  simplifications  are  necessary.  ...  [One]  approach, 

...,  is  to  use  parallel  machines,  each  handling  only  a  part  of  the  sim¬ 
ulated  system  ...[Another]  approach,  ...,  is  to  move  to  a  higher  and 
much  faster  level  of  modeling,  the  event-driven  or  “state-machine,” 
scheme.  In  order  to  retain  the  essential  features  of  the  modeled  sys¬ 
tem,  the  construction  of  this  level  of  representation  should  be  based 
heavily  on  the  functional  rules  that  were  formulated  as  a  result  of 
detailed  exploration  of  the  [detailed]  models. 

Aside  from  the  advantageous  information  yield  of  working  with  biophysical 
models,  we  were  attracted  to  this  complementary  modeling  paradigm  because 
much  of  the  detailed  data  needed  for  significant  aspects  of  the  ventricular  sim¬ 
ulation  study  was  unavailable  from  the  research  hterature  and  the  (unproven) 
possibility  that  reasonable  hypotheses  about  the  missing  data  could  be  gen¬ 
erated  through  biophysical  simulations.  Among  our  projected  areas  of  study 
for  this  purpose  were  the  action  potential  under  myocardial  ischemia  and  the 
self-organization  of  refractory  period  distributions  and  “cardiac  memory”  under 
various  scenarios  of  conduction  pathology  (such  as  bundle  blocks).  However,  the 
extent  of  our  work  in  this  area  reached  only  so  far  as  simulation  of  fairly  simple 
ionic  models  on  small  two-dimensional  arrays. 

We  have  tried  to  describe  our  view  of  a  proper  framework  for  modeling  and 
simulation  of  bioelectric  phenomena;  next,  we  would  like  to  describe  a  formal 
framework  for  model-based  derivation  of  diagnostic  criteria. 

4.2  The  Work  Unit  in  a  Broader  Context 

As  part  of  our  efforts  to  understand  how  we  could  best  enhance  the  ECG,  we 
looked  at  the  large-scale  structure  of  the  study  design.  The  study  could  roughly 
be  broken  down  into  two  components:  the  model-based  generation  of  hypotheses 
about  optimal  lead  systems  and  diagnostic  criteria,  and  clinical  investigation 
and  validation  of  these  hypotheses.  While  the  clinical  aspects  of  the  study  were 
clearly  outside  the  scope  of  our  in-house  efforts,  it  became  our  view  that  the 
model-based  hypothesis  generation  component  of  the  study  could  have  benefited 
from  a  more  formal  approach  to  its  design. 

In  particular,  ECG  diagnostic  criteria  of  the  type  sought  in  the  study  ba¬ 
sically  consist  of  rules  for  interpreting  the  electric  potential  measurements  in 
various  combinations  of  ECG  leads  (e.g.,  see  [16]).  Ischemic  electrophysiology 
can  be  conceptually  modeled  as  a  perturbation  of  normal  electrophysiology;  di¬ 
agnostic  criteria  for  ischemia  can  be  viewed  as  an  attempt  to  create  a  mapping 
between  perturbations  from  normal  in  ECG  potentials  and  perturbations  from 
normal  in  the  underlying  electrophysiology.  The  components  of  this  mapping 
provide  the  proper  abstractions  to  organize  supporting  modeling  and  simulation 
studies.  This  mapping  can  be  formulated  as  a  set  of  pairings  between  combi- 
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nations  of  perturbations  in  the  electrophysiology  of  various  coronary  perfusion 
fields  and  responses  in  ECG  lead  sets  or  feature  sets  of  interest. 

Selvester  and  Solomon’s  origincJ  study  design  to  enhance  ECG  diagnostics 
for  ischemia  included  the  generation  of  simulated  ECG  data  for  “all  possible 
combinations”  of  ischemia  at  different  levels  of  severity  in  all  of  the  coronary 
perfusion  fields  [1].  This  was  to  be  followed  by  an  analysis  of  the  resulting 
simulation  database  to  determine  optimal  electrode  sites  and  diagnostic  criteria 
for  various  combinations  of  ischemic  perfusion  fields. 

Beside  arguments  concerning  the  scientific  feasibility  of  specifying  ischemic 
effects  in  the  detail  required  by  Selvester  and  Solomon’s  study  design,  we  could 
see  potential  methodological  problems.  In  particular,  there  was  no  specified 
approach  for  analyzing  the  simulation  database,  nor  was  there  any  formalism 
for  managing  the  generation  of  simulation  data  in  terms  of  the  desired  framework 
of  abstractions.  A  natural  approach  to  these  types  of  problems  can  be  found 
in  the  field  of  model-based  diagnosis  (a  branch  of  the  field  of  expert  systems 
technology).  The  advantage  of  making  this  conceptual  association  is  that  the 
methods  of  model-based  diagnosis  are  usually  formulated  directly  in  terms  of 
the  abstractions  found  at  the  diagnosis  level.  The  models  referred  to  in  the  term 
“model-based”  are  usually  qualitative  models  describing  the  causal  behavior  of 
the  system  under  diagnosis  in  terms  of  appropriate  diagnostic  abstractions. 

An  excellent  prototype  of  this  approach  is  found  in  the  work  performed 
by  Bratko,  et  al.,  on  their  KARDIO  system  for  automated  diagnosis  of  cardiac 
arrhythmias  [18].  Their  work  involved  comprehensive  simulation  of  many  classes 
of  arrhythmias,  resulting  in  a  large  database  of  pairings  between  conditions  in 
their  heart  model  and  corresponding  ECG  features  in  the  form  of  logic  formulas. 
Formal  machine  learning  techniques  were  then  applied  to  the  database  to  extract 
diagnostic  rules. 

We  believe  that  the  simulation  of  basic  electrophysiological  scenarios  can  be 
combined  with  perturbative  sensitivity  analysis  to  derive  a  qualitative  causal 
model  of  the  ECG  for  ischemic  heart  conditions.  The  qualitative  model  would 
then  provide  a  translation  between  fundamental  mathematical  models  and  the 
types  of  symbolic  abstractions  which  can  be  utdized  by  model-based  diagnosis 
techidques.  This  would  facilitate  the  use  of  formal  machine  learning  techniques 
in  analyzing  such  simulation  results  and  relating  them  in  a  weU-defined  way  to 
diagnostic  criteria. 

The  work  of  Selvester  and  Solomon  on  criteria  for  infarction  and  ischemia 
has  shown  that  simulations  studies  can  be  used  to  help  understand  the  structure 
of  information  in  the  ECG  in  terms  of  underlying  biophysical  principles  and  how 
this  information  can  be  used  to  enhance  the  diagnostic  capabilities  of  the  ECG. 
The  ECG  enhancement  project  as  a  whole,  including  both  the  contracted  effort 
and  this  work  unit,  can  be  viewed,  within  the  context  of  this  emerging  discipline 
of  model-based  enhancement  of  diagnostic  systems,  as  an  informal  approach  to 
“model-based  knowledge  acquisition”  [18]. 

The  review  of  Uckun  [19]  gives  a  broader  overview  of  applications  of  model- 
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based  reasoning  in  biomedicine,  and  the  paper  of  Tong  and  Widman  [20]  pro¬ 
vides  a  recent  example  of  the  application  of  this  paradigm  to  ECG  diagnostics. 
The  paper  of  Siregar,  et  al.,  [21]  indicates  hovir  more  general  models  can  fit  into 
this  framework,  as  well  as  a  concise  and  lucid  review  of  model-based  diagno¬ 
sis.  As  the  literature  shows,  model-based  reasoning  is  currently  a  very  active 
subfield  of  artificial  intelligence  research,  with  potential  for  application  to  any 
diagnostic  domain  in  which  the  system  under  diagnosis  can  be  modeled  by  either 
quantitative  or  qualitative  methods. 

We  also  want  to  point  out  that  a  contemporary  approach  to  the  type  of 
simulation  study  undertaken  in  the  ECG  enhancement  effort  would  probably 
be  quite  different  than  that  initiated  in  1984.  The  complexities  arising  from  the 
spatial  distribution  of  parameters  and  the  irregular  geometries  of  anatomical 
domains  would  be  addressed  through  the  application  of  high-performance  com¬ 
puters  and  tools  such  as  magnetic  resonance  imaging  (MRI),  computer-aided 
design,  and  the  finite  element  method  ([22],  [23]).  The  use  of  such  tools  is  cur¬ 
rently  revolutionizing  the  discipline  of  biomediccJ  modeling  and  simulation  [22]. 
The  recent  work  of  MUler,  et  al.,  [24]  may  hint  at  a  future  where  the  anatomy 
and  physiology  of  populations  of  individuals  can  be  computationally  modeled, 
and  analyzed  by  machine-oriented  techniques,  to  achieve  an  understanding  of 
complexities  that  humans  could  not  otherwise  hope  to  comprehend. 

4.3  Lessons  Learned 

In  conclusion,  we  mention  a  few  lessons  learned  from  our  general  experience  on 
this  project. 

Computationrd  methods  for  the  simulation  and  analysis  of  spatially 
complex  biological  systems.  There  are  significant  computational  difficul¬ 
ties  involved  in  simulating  complex  physiological  interactions  in  systems  where 
spaticil  extent  is  large  vdth  respect  to  the  scale  at  which  the  biophysical  in¬ 
teractions  take  place.  Such  problems  can  be  computationally  intractable  un¬ 
less  simplifications  can  be  made  in  the  simulation  requirements.  Methods  are 
needed  which  optimally  preserve  details  of  the  underlying  model  versus  compu¬ 
tational  requirements.  Such  methods  could  range  over  many  different  combi¬ 
nations  of  simulation  software  technique  and  high-performance  computing.  In 
dealing  with  systems  whose  components  possess  inhomogeneous  characteristics, 
and  are  changing  in  various  scenarios  of  interest,  it  is  important  to  use  methods 
which  are  not  only  computationally  efficient,  but  also  flexible.  Finally,  when  the 
input  set  to  a  simrdation  includes  complex  spatial  distributions  of  parameters, 
and  the  output  set  is  of  a  similar  nature,  preparation  and  analysis  of  a  system¬ 
atic  simulation  study  can  become  extremely  cumbersome.  Thus,  there  is  a  great 
need  for  computational  tools  allowing  easy  interactive  generation,  manipulation, 
and  analysis  of  data  defined  over  geometric  domains  in  three  dimensions. 


19 


The  uniqueness  of  Air  Force  requirements.  In  recent  years,  there  has 
been  considerable  activity  in  the  modeling  and  simulation  of  cardiac  electro¬ 
physiology  by  the  mathematical  biology  and  biomedical  engineering  commu¬ 
nities.  Most  of  this  research,  however,  has  been  directed  toward  modeling  the 
propagation  of  the  depolarization  wavefront  in  simplified  models  for  the  purpose 
of  gaining  insight  into  the  nature  of  intra-myo cardial  arrhythmias. 

Our  interests,  in  the  EGG  enhancement  effort,  included  those  of  the  broader 
community,  but  extended  beyond  them.  We  were  primarily  interested  in  gaining 
insight  into  the  detection  of  asymptomatic  coronary  artery  disease  via  electro¬ 
cardiography,  as  well  as  related  issues  such  as  modeling  the  electro-mechaiucal 
interaction.  These  objectives  direct  our  attention  to  what  is  happening  behind 
the  wavefront  of  myocardial  depolarization,  for  it  is  there  that  the  most  sensitive 
indicators  of  asymptomatic  CAD  arise,  and  it  is  there  that  the  fine  details  of 
electro-mechanical  interaction  occur.  Moreover,  unlike  much  of  the  work  on  the 
modeling  of  arrhythmias,  our  program  objectives  required  that  we  address  the 
modehng  of  both  the  depolarization  and  repolarization  processes  over  the  entire 
muscle  mass  of  the  ventricles,  with  attention  to  the  details  of  physiologically 
relevant  inhomogeneities  as  they  change  in  various  scenarios  of  interest,  such  as 
rest,  exercise,  and  -f-Gz  stress. 

These  areas  of  emphasis  stand  apart  from  the  concerns  of  the  EGG  and  elec¬ 
trophysiology  modeling  communities  at-large,  and  constitutes  another  demon¬ 
stration  that  the  Air  Force  presents  requirements  for  biomedical  investigation 
which  may  not  be  directly  addressed  by  the  general  scientific  community. 
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